setwd("c:/YourPath/YourFolder")

if (!require("pacman")) install.packages("pacman")
pacman::p_load(tidyverse, readstata13, dplyr, ggplot2, lubridate, stringr, scales, ggrepel, lubridate, xlsx,
               MASS, mgcv, tidygam, flextable, modelsummary, foreign, marginaleffects, haven, car, ggpubr,
               gridExtra)

#read data
emendamenti_sub <- read.csv("emendamenti_sub.csv", header = TRUE)
emendamenti_sub <- emendamenti_sub[,-1]
emendamenti_sub$opposition<- as.factor(emendamenti_sub$opposition)


emendamenti_sub2 <- subset(emendamenti_sub, fiducia==0)
emendamenti_sub3 <- subset(emendamenti_sub, camera_pl==1)
emendamenti_sub4 <- subset(emendamenti_sub3, fiducia==0)

# model
##########################################
# variables' name:                       #
# opposition1 = opposition               #
# diff1 = Ideological distance           #
# size_rel = Seat share                  #
# gov_exp = Prior exp. in gov.           #
# populist = Populist party              #
# gov_dist = Gov. policy divisions       #
# camera_pl = First reading Chamber      #
# commissioni = More than one committee  #
# prop_Msalute = Ministry of Health      #
##########################################

m1_qp1 <- gam(em_tot ~ opposition + diff1 + size_rel + gov_exp  + populist +
                + gov_dist + 
                + camera_pl +  commissioni + prop_Msalute +
                s(days, by=opposition, fx = FALSE, k=10), 
              data = emendamenti_sub, family = quasipoisson, method = "REML")

m1_qp2 <- gam(em_tot ~ opposition + diff1 + size_rel + gov_exp  + populist + 
                + gov_dist +
                + camera_pl + commissioni + prop_Msalute +
                s(days, by=opposition, fx = FALSE, k=10), 
              data = emendamenti_sub2, family = quasipoisson, method = "REML")

m1_qp3 <- gam(em_tot ~ opposition + diff1 + size_rel + gov_exp  + populist +
                + gov_dist +
                +  commissioni + prop_Msalute +
                s(days, by=opposition, fx = FALSE, k=10), 
              data = emendamenti_sub3, family = quasipoisson, method = "REML")

m1_qp4 <- gam(em_tot ~ opposition + diff1 + size_rel + gov_exp  + populist +
                + gov_dist +
                +  commissioni + prop_Msalute +
                s(days, by=opposition, fx = FALSE, k=10), 
              data = emendamenti_sub4, family = quasipoisson, method = "REML")


# Table 1
m1 <- flextable::as_flextable(m1_qp1)
m2 <- flextable::as_flextable(m1_qp2)
m3 <- flextable::as_flextable(m1_qp3)
m4 <- flextable::as_flextable(m1_qp4)


# Figure 1
fig1data <- data.frame(
  mese = c("Sep-19", "Oct-19", "Dec-19", "Feb-20", "Mar-20", "Apr-20", "Jun-20", "Aug-20", "Oct-20", "Dec-20"),
  approvazione = c(44, 40, 42, 44, 71, 63, 60, 57, 55, 57)
)

fig1data$mese <- factor(fig1data$mese, 
                        levels = c("Sep-19", "Oct-19", "Dec-19", "Feb-20", "Mar-20", "Apr-20", "Jun-20", "Aug-20", "Oct-20", "Dec-20"))

ggplot(fig1data, aes(x = mese, y = approvazione)) +
  geom_bar(stat = "identity") +
  theme_minimal() +
  labs(title = "",
       x = "Month/Year",
       y = "Approval rating (%)") +
  scale_x_discrete(guide = guide_axis(angle = 90)) +
  theme(text = element_text(size=16))

# Figure 2
emendamenti_sub_sum <- emendamenti_sub %>% 
  group_by(months) %>% summarise(tutti =n_distinct(ac), fiducia = n_distinct(fiducia))

sum(emendamenti_sub_sum$tutti)
sum(emendamenti_sub_sum$fiducia)

emendamenti_sub_sum$months <- as.factor(emendamenti_sub_sum$months) 
emendamenti_sub_sum$netto <- emendamenti_sub_sum$tutti - emendamenti_sub_sum$fiducia

emendamenti_sub_sum$tutti <- NULL

emendamenti_sub_sum_long <- gather(emendamenti_sub_sum, key = "Type", value = "sum", -months)

ggplot(emendamenti_sub_sum_long, aes(x=months, y=sum, fill=Type)) + 
  geom_bar(stat = "identity") +
  theme_minimal() +
  labs(title="", x="Year/Month", y="N. of Decree-Law Coversion Bills") +
  scale_x_discrete(guide = guide_axis(angle = 90)) +
  guides(fill = guide_legend(title = "")) +
  theme(legend.position = "bottom") +
  scale_fill_manual(values = c("#999999", "#000000"), labels = c("With Confidence", "Without Confidence")) +
  theme(text = element_text(size=16))

# Figure 3
p1 <- plot_slopes(m1_qp1, variables = 'opposition',
                  condition = c('days'),
                  type = 'link', rug=TRUE) +
  annotate("rect", xmin=709, xmax=754, ymin = -Inf, ymax = Inf,
           alpha = .4) +
  annotate("text", x = 730, y = 4, label = "1st \n wave") +
  annotate("rect", xmin=937, xmax=998, ymin = -Inf, ymax = Inf,
           alpha = .4) +
  annotate("text", x = 960, y = 6, label = "2nd \n wave") +
  annotate("rect", xmin=1074, xmax=1119, ymin = -Inf, ymax = Inf,
           alpha = .4) +
  annotate("text", x = 1090, y = 8, label = "3rd \n wave") +
  geom_hline(yintercept=0, linetype="dashed", 
             color = "black", size=0.5) +
  ggtitle("M1") +
  scale_y_continuous(limits = c(0, 10))

p2 <- plot_slopes(m1_qp2, variables = 'opposition',
                  condition = c('days'),
                  type = 'link', rug=TRUE) +
  annotate("rect", xmin=709, xmax=754, ymin = -Inf, ymax = Inf,
           alpha = .4) +
  annotate("text", x = 730, y = 4, label = "1st \n wave") +
  annotate("rect", xmin=937, xmax=998, ymin = -Inf, ymax = Inf,
           alpha = .4) +
  annotate("text", x = 960, y = 6, label = "2nd \n wave") +
  annotate("rect", xmin=1074, xmax=1119, ymin = -Inf, ymax = Inf,
           alpha = .4) +
  annotate("text", x = 1090, y = 8, label = "3rd \n wave") +
  geom_hline(yintercept=0, linetype="dashed", 
             color = "black", size=0.5) +
  ggtitle("M2") +
  scale_y_continuous(limits = c(0, 10))

p3 <- plot_slopes(m1_qp3, variables = 'opposition',
                  condition = c('days'),
                  type = 'link', rug=TRUE) +
  annotate("rect", xmin=709, xmax=754, ymin = -Inf, ymax = Inf,
           alpha = .4) +
  annotate("text", x = 730, y = 2, label = "1st \n wave") +
  annotate("rect", xmin=937, xmax=998, ymin = -Inf, ymax = Inf,
           alpha = .4) +
  annotate("text", x = 960, y = 4, label = "2nd \n wave") +
  annotate("rect", xmin=1074, xmax=1119, ymin = -Inf, ymax = Inf,
           alpha = .4) +
  annotate("text", x = 1090, y = 6, label = "3rd \n wave") +
  geom_hline(yintercept=0, linetype="dashed", 
             color = "black", size=0.5) +
  ggtitle("M3") +
  scale_y_continuous(limits = c(0, 10))

p4 <- plot_slopes(m1_qp4, variables = 'opposition',
                  condition = c('days'),
                  type = 'link', rug=TRUE) +
  annotate("rect", xmin=709, xmax=754, ymin = -Inf, ymax = Inf,
           alpha = .4) +
  annotate("text", x = 730, y = 3, label = "1st \n wave") +
  annotate("rect", xmin=937, xmax=998, ymin = -Inf, ymax = Inf,
           alpha = .4) +
  annotate("text", x = 960, y = 5, label = "2nd \n wave") +
  annotate("rect", xmin=1074, xmax=1119, ymin = -Inf, ymax = Inf,
           alpha = .4) +
  annotate("text", x = 1090, y = 7, label = "3rd \n wave") +
  geom_hline(yintercept=0, linetype="dashed", 
             color = "black", size=0.5) +
  ggtitle("M4") +
  scale_y_continuous(limits = c(0, 10))

ggarrange(p1, p2, p3, p4, ncol = 2, nrow = 2) 

# Figure 4
p1c <- plot_comparisons(m1_qp1, variables = 'opposition',
                        condition = c('days'),
                        type = 'response', rug=TRUE) +
  annotate("rect", xmin=709, xmax=754, ymin = -Inf, ymax = Inf,
           alpha = .4) +
  annotate("text", x = 730, y = 15, label = "1st \n wave") +
  annotate("rect", xmin=937, xmax=998, ymin = -Inf, ymax = Inf,
           alpha = .4) +
  annotate("text", x = 960, y = 20, label = "2nd \n wave") +
  annotate("rect", xmin=1074, xmax=1119, ymin = -Inf, ymax = Inf,
           alpha = .4) +
  annotate("text", x = 1090, y = 25, label = "3rd \n wave") +
  geom_hline(yintercept=0, linetype="dashed", 
             color = "black", size=0.5) +
  ggtitle("M1") +
  scale_y_continuous(limits = c(0, 100))

p2c <- plot_comparisons(m1_qp2, variables = 'opposition',
                        condition = c('days'),
                        type = 'response', rug=TRUE) +
  annotate("rect", xmin=709, xmax=754, ymin = -Inf, ymax = Inf,
           alpha = .4) +
  annotate("text", x = 730, y = 15, label = "1st \n wave") +
  annotate("rect", xmin=937, xmax=998, ymin = -Inf, ymax = Inf,
           alpha = .4) +
  annotate("text", x = 960, y = 20, label = "2nd \n wave") +
  annotate("rect", xmin=1074, xmax=1119, ymin = -Inf, ymax = Inf,
           alpha = .4) +
  annotate("text", x = 1090, y = 25, label = "3rd \n wave") +
  geom_hline(yintercept=0, linetype="dashed", 
             color = "black", size=0.5) +
  ggtitle("M2") +
  scale_y_continuous(limits = c(0, 100))

p3c <- plot_comparisons(m1_qp3, variables = 'opposition',
                        condition = c('days'),
                        type = 'response', rug=TRUE) +
  annotate("rect", xmin=709, xmax=754, ymin = -Inf, ymax = Inf,
           alpha = .4) +
  annotate("text", x = 730, y = 15, label = "1st \n wave") +
  annotate("rect", xmin=937, xmax=998, ymin = -Inf, ymax = Inf,
           alpha = .4) +
  annotate("text", x = 960, y = 20, label = "2nd \n wave") +
  annotate("rect", xmin=1074, xmax=1119, ymin = -Inf, ymax = Inf,
           alpha = .4) +
  annotate("text", x = 1090, y = 25, label = "3rd \n wave") +
  geom_hline(yintercept=0, linetype="dashed", 
             color = "black", size=0.5) +
  ggtitle("M3") +
  scale_y_continuous(limits = c(-1, 100))

p4c <- plot_comparisons(m1_qp4, variables = 'opposition',
                        condition = c('days'),
                        type = 'response', rug=TRUE) +
  annotate("rect", xmin=709, xmax=754, ymin = -Inf, ymax = Inf,
           alpha = .4) +
  annotate("text", x = 730, y = 15, label = "1st \n wave") +
  annotate("rect", xmin=937, xmax=998, ymin = -Inf, ymax = Inf,
           alpha = .4) +
  annotate("text", x = 960, y = 20, label = "2nd \n wave") +
  annotate("rect", xmin=1074, xmax=1119, ymin = -Inf, ymax = Inf,
           alpha = .4) +
  annotate("text", x = 1090, y = 25, label = "3rd \n wave") +
  geom_hline(yintercept=0, linetype="dashed", 
             color = "black", size=0.5) +
  ggtitle("M4") +
  scale_y_continuous(limits = c(0, 100))

ggarrange(p1c, p2c, p3c, p4c, ncol = 2, nrow = 2)

# Appendix

# Table A1
party <- emendamenti_sub %>% 
  group_by(partito) %>%
  dplyr::summarize(mean = mean(em_tot, na.rm=TRUE),
                   max = max(em_tot, na.rm=TRUE),
                   SD = sd(em_tot, na.rm=TRUE),
                   n = n())

# Table A2
emendamenti_sub$populist <- as.factor(emendamenti_sub$populist)

pp = gam(em_tot ~  size_rel + gov_exp  + 
           + gov_dist + 
           + camera_pl +  prop_Msalute + commissioni + populist +
           s(days, by=diff1, fx = FALSE, k=10) +
           s(days, by=populist, fx = FALSE, k=10), 
         data = emendamenti_sub, family = quasipoisson, method = "REML")

a2 <- flextable::as_flextable(pp)

# Figure A1
plot_slopes(pp, variables = 'diff1',
            condition = c('days'),
            type = 'link', rug=TRUE) +
  annotate("rect", xmin=709, xmax=754, ymin = -Inf, ymax = Inf,
           alpha = .4) +
  annotate("text", x = 730, y = .2, label = "1st \n wave") +
  annotate("rect", xmin=937, xmax=998, ymin = -Inf, ymax = Inf,
           alpha = .4) +
  annotate("text", x = 960, y = .4, label = "2nd \n wave") +
  annotate("rect", xmin=1074, xmax=1119, ymin = -Inf, ymax = Inf,
           alpha = .4) +
  annotate("text", x = 1090, y = .6, label = "3rd \n wave") +
  geom_hline(yintercept=0, linetype="dashed", 
             color = "black", size=0.5) 

# Figure A2
plot_slopes(pp, variables = 'populist',
            condition = c('days'),
            type = 'link', rug=TRUE) +
  annotate("rect", xmin=709, xmax=754, ymin = -Inf, ymax = Inf,
           alpha = .4) +
  annotate("text", x = 730, y = .2, label = "1st \n wave") +
  annotate("rect", xmin=937, xmax=998, ymin = -Inf, ymax = Inf,
           alpha = .4) +
  annotate("text", x = 960, y = .4, label = "2nd \n wave") +
  annotate("rect", xmin=1074, xmax=1119, ymin = -Inf, ymax = Inf,
           alpha = .4) +
  annotate("text", x = 1090, y = .6, label = "3rd \n wave") +
  geom_hline(yintercept=0, linetype="dashed", 
             color = "black", size=0.5) 

